import os
import sys


#-- general import
import os
#-- OGR 
try:
    from osgeo import ogr
except ImportError:
    import ogr

#-- Shapely
from shapely.wkb import loads as wkbloads
from shapely.wkt import loads as wktloads
from shapely.geometry import Point, LineString, Polygon
from shapely.ops import cascaded_union
'''
Created on 12/01/2011

@author: jguerrer
'''
def writeBackShapeFile(eleFile, nodeFile):
    # the elefile contains the actual triangles
    # the nodefile contains the list of vertexes or points
    
    #first create a dictionary of nodefiles
    
    
    print "Reading node file: " + nodeFile
    f=open(nodeFile,'r')
    headerNode = f.readline()
    elements = headerNode.split()
    
    nodesNum = elements[0];
    nodesDim = elements[1];
    #skip the next elements like the attributes

    nodes = {}
    data=f.readlines()
    for x in xrange(len(data)-1):#for each line
        d = data[x]
        #print d
        elements = d.split()
        # pt_number  x_coordinate  y_coordinate z_attribute
        #print "0:" + elements[0] + " 1:" + str(elements[1]) + " 2:" + str(elements[2]) + " 3:" + str(elements[3])
        
        nodes[elements[0]] = { 'id':elements[0] ,'x':elements[1] , 'y':elements[2], 'z':elements[3] }
        
    #now points are ready to use    
    
    print "Reading ele file: " + eleFile
    f=open(eleFile,'r')
    headerNode = f.readline()
    elements = headerNode.split()
    triangles = [];
    data=f.readlines()
    for x in xrange(len(data)-1):#for each line
        d = data[x]
        #print d
        elements = d.split()
        p1 = elements[1]
        p2 = elements[2]
        p3 = elements[3]
     
        
       ############ triangles
        wkt = "POLYGON(( "   + \
        nodes[p1]['x'] + " " + nodes[p1]['y']  + " , "  + \
        nodes[p2]['x'] + " " + nodes[p2]['y']  + " , "  + \
        nodes[p3]['x'] + " " + nodes[p3]['y']  + " , "  + \
        nodes[p1]['x'] + " " + nodes[p1]['y']  + "))";
        #print wkt;
        triangles.append({'geom':wkt });
    ##start writing a shapefile os send it to a method
    return triangles;
def write_shapefile(lsFeat, outfile):
    """
    Create a new shapefile and writes it to the disk.
    
    Input:
        lsFeat:    list of Features (a dico with geom/name/pop)
        outfile:   path for the new file to create
    Output:
        - (the shapefile is written to the disk)

    """
    driver = ogr.GetDriverByName('ESRI Shapefile') #-- we create a new SHP file
    if os.path.exists(outfile):
            driver.DeleteDataSource(outfile) #-- if it exists, overwrite it
    ds = driver.CreateDataSource(outfile)
    layer = ds.CreateLayer(outfile, geom_type=ogr.wkbPolygon) #-- we create a SHP with polygons
    
    #-- create the 1 attributes (called Fields in OGR jargon) "name" and "population"
    fd = ogr.FieldDefn()
    fd.SetName('height')
    fd.SetType(ogr.OFTString)  #-- type is String
    layer.CreateField(fd)
        
    for i in lsFeat:
        f = ogr.Feature(feature_def=layer.GetLayerDefn())
        
        p = ogr.CreateGeometryFromWkt(i['geom'])
        f.SetGeometry(p)
#        f.SetField(0, i['height'])
        f.SetField(0, 0)
        layer.CreateFeature(f)
        f.Destroy()
    ds.Destroy()
    print "\nShapefile saved as:", outfile 
    
def main():
    nodefile = 'test.1.node'
    elefile = 'test.1.ele'
    outputShapefile = 'outputTriangulation.shp'
    triangles = writeBackShapeFile(elefile,nodefile);
    write_shapefile(triangles, outputShapefile)
    
    pass;

if __name__ == '__main__':
    main();

